#############################
#####MULTILEVEL ANALYSIS#####
#############################

#MARIANA CARMO DUARTE

#AUGUST 2024

#Please feel free to contact me at mariana.carmoduarte@ics.ulisboa.pt in case
#you have any doubts about the code and data. 

#GENERAL LIBRARIES
library(foreign)
library(tidyverse)
library(dplyr)
library(ggplot2)
library(viridis)
library(viridisLite)
library(dplyr)
library(fastDummies)
library(tidyverse)
library(ggplot2)
library(MASS)
library(calibrate)
library(psych)
library(lattice)
library(car)
library(plyr)
library(gplots)
library(hrbrthemes)
library(rms)
library(tidyverse)
library(MASS)
library(calibrate)
library(psych)
library(lattice)
library(plyr)
library(dplyr)
library(stargazer)
library(xtable)
library(foreign)
library(texreg)
library(readstata13)
library(ggplot2)
library(tidyr)
library(moments)
library(lmtest)
library(estimatr)
library(olsrr)
library(interplot)
library(sjmisc)
library(sjPlot)
library(interplot)
library(gmodels)
library(oddsratio)
library(epiR)
library(Rmisc)
library(Zelig)
library(mfx)
library(car)
library(skimr)
library(prediction)
library(sandwich)
library(gridExtra)
library(irr)
library(miceadds)
library(stargazer)
library(clusterSEs)
library(lmtest)
library(texreg)
library(margins)
library(coefplot)

#DATA MANIPULATION
library(haven)       
library(dplyr)       
library(psych)       

#MODELING
library(lme4)        
library(nlme)        
library(arm)         
library(optimx)      
library(performance)  
library(insight)     
library(ordinal)
library(interflex)

#VISUALIZATION
library(ggplot2)     
library(scales)      
library(interplot)   
library(sjPlot)      
library(rockchalk)   
library(sjlabelled)  
library(ggeffects)   

###
library(HLMdiag)
library(DHARMa)
library(jtools)

#LOAD INDIVIDUAL-LEVEL ADATA
setwd("")
ESS_all <- read.spss("DATA_FINAL.sav", use.value.labels = TRUE, to.data.frame = TRUE)

#PREPARATION OF VARIABLES AND CENTERING
#NA's and NUMERIC
#INDEPENDENT VARIABLES
ESS_all$immig_index <- as.numeric(ESS_all$immig_index)
ESS_all$trstn_index <- as.numeric(ESS_all$trstn_index)
ESS_all$imptrad <- as.numeric(ESS_all$imptrad)
ESS_all$ipstrgv <- as.numeric(ESS_all$ipstrgv)
ESS_all$polarization_seats.x <- as.numeric(ESS_all$polarization_seats.x)
ESS_all$eduyrs <- as.numeric(ESS_all$eduyrs)
ESS_all$edulvlb <- as.numeric(ESS_all$edulvlb)
ESS_all$hinctnta <- as.numeric(ESS_all$hinctnta)
ESS_all$hincfel <- as.numeric(ESS_all$hincfel)
ESS_all$agea <- as.numeric(ESS_all$agea)
ESS_all$stfeco <- as.numeric(ESS_all$stfeco)
ESS_all$stfgov <- as.numeric(ESS_all$stfgov)
ESS_all$trstprl <- as.numeric(ESS_all$trstprl)
ESS_all$pol_extrem <- as.numeric(ESS_all$pol_extrem)
ESS_all$hinctnta <- as.numeric(ESS_all$hinctnta)
ESS_all$hincfel <- as.numeric(ESS_all$hincfel)
ESS_all$imdfetn <- as.numeric(ESS_all$imdfetn)
ESS_all$lrscale <- as.numeric(ESS_all$lrscale)
ESS_all$stfdem <- as.numeric(ESS_all$stfdem)
ESS_all$eisced <- as.numeric(ESS_all$eisced)
ESS_all$eu_balance <- as.numeric(ESS_all$eu_balance)
ESS_all$edulvla <- as.numeric(ESS_all$edulvla)
ESS_all$imig_blive <- as.numeric(ESS_all$imig_blive)
ESS_all$yearsince <- as.numeric(ESS_all$yearsince)
ESS_all$gincdif <- as.numeric(ESS_all$gincdif)
ESS_all$rlgdgr <- as.numeric(ESS_all$rlgdgr)
ESS_all$euftf <- as.numeric(ESS_all$euftf)
ESS_all$trstn_index <- as.numeric(ESS_all$trstn_index)
ESS_all$rural <- as.numeric(ESS_all$rural)
ESS_all$male <- as.numeric(ESS_all$male)


#CENTERING
centFUN <- function(x) {
  x - mean(x, na.rm = TRUE)
}

#GRAND-MEAN CENTERING
ESS_all$immig_indexCGM <- centFUN(ESS_all$immig_index)
ESS_all$trstn_indexCGM <- centFUN(ESS_all$trstn_index)
ESS_all$eduyrsCGM <- centFUN(ESS_all$eduyrs)
ESS_all$edulvlbCGM <- centFUN(ESS_all$edulvlb)
ESS_all$hinctntaCGM <- centFUN(ESS_all$hinctnta)
ESS_all$hincfelCGM <- centFUN(ESS_all$hincfel)
ESS_all$ageaCGM <- centFUN(ESS_all$agea)
ESS_all$stfecoCGM <- centFUN(ESS_all$stfeco)
ESS_all$stfgovCGM <- centFUN(ESS_all$stfgov)
ESS_all$trstprlCGM <- centFUN(ESS_all$trstprl)
ESS_all$pol_extremCGM <- centFUN(ESS_all$pol_extrem)
ESS_all$hinctntaCGM <- centFUN(ESS_all$hinctnta)
ESS_all$hincfelCGM <- centFUN(ESS_all$hincfel)
ESS_all$imdfetnCGM <- centFUN(ESS_all$imdfetn)
ESS_all$lrscaleCGM <- centFUN(ESS_all$lrscale)
ESS_all$stfdemCGM <- centFUN(ESS_all$stfdem)
ESS_all$eiscedCGM <- centFUN(ESS_all$eisced)
ESS_all$eu_balanceCGM <- centFUN(ESS_all$eu_balance)
ESS_all$edulvlaCGM <- centFUN(ESS_all$edulvla)
ESS_all$imig_bliveCGM <- centFUN(ESS_all$imig_blive)
ESS_all$trstepCGM <- centFUN(ESS_all$trstep)
ESS_all$gincdifCGM <- centFUN(ESS_all$gincdif)
ESS_all$rlgdgrCGM <- centFUN(ESS_all$rlgdgr)
ESS_all$ruralCGM <- centFUN(ESS_all$rural)
ESS_all$maleCGM <- centFUN(ESS_all$male)

#GROUP-MEAN CENTERING
ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(trstepCWC = centFUN(trstep))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(gincdifCWC = centFUN(gincdif))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(pol_extremCWC = centFUN(pol_extrem))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(imig_bliveCWC = centFUN(imig_blive))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(lrscaleCWC = centFUN(lrscale))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(edulvlaCWC = centFUN(edulvla))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(eiscedCWC = centFUN(eisced))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(immig_indexCWC = centFUN(immig_index))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(imdfetnCWC = centFUN(imdfetn))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(trstn_indexCWC = centFUN(trstn_index))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(trstprlCWC = centFUN(trstprl))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(eduyrsCWC = centFUN(eduyrs))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(ruralCWC = centFUN(rural))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(maleCWC = centFUN(male))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(edulvlbCWC = centFUN(edulvlb))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(hinctntaCWC = centFUN(hinctnta))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(hincfelCWC = centFUN(hincfel))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(ageaCWC = centFUN(agea))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(stfecoCWC = centFUN(stfeco))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(stfgovCWC = centFUN(stfgov))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(stfdemCWC = centFUN(stfdem))

ESS_all <- ESS_all %>%
  group_by(election_id) %>%
  mutate(rlgdgrCWC = centFUN(rlgdgr))


#Grand-mean L2 Predictors
#L2 variable value in each L2 group

temp_ess<- ESS_all%>%
  group_by(election_id) %>%
  summarise(politicisation.x = mean(politicisation.x, na.rm = TRUE))

temp_esss<- ESS_all%>%
  group_by(election_id) %>%
  summarise(polarization_seats.x = mean(polarization_seats.x, na.rm = TRUE))

temp_ess1<- ESS_all%>%
  group_by(election_id) %>%
  summarise(eu_polarisation_ps.x = mean(eu_polarisation_ps.x, na.rm = TRUE))

temp_esss2<- ESS_all%>%
  group_by(election_id) %>%
  summarise(eu_salience_ps.x = mean(eu_salience_ps.x, na.rm = TRUE))

temp_essy<- ESS_all%>%
  group_by(election_id) %>%
  summarise(politicisation.y = mean(politicisation.y, na.rm = TRUE))

temp_ess1y<- ESS_all%>%
  group_by(election_id) %>%
  summarise(eu_polarisation_ps.y = mean(eu_polarisation_ps.y, na.rm = TRUE))

temp_esss2y<- ESS_all%>%
  group_by(election_id) %>%
  summarise(eu_salience_ps.y = mean(eu_salience_ps.y, na.rm = TRUE))

temp_essssss<- ESS_all%>%
  group_by(election_id) %>%
  summarise(yearsince = mean(yearsince, na.rm = TRUE))

#subtract the mean of the L2 variable value across country-year from
# the country-year predictor values
#Compute the mean of POLITICISATION
mean(temp_ess$politicisation.x, na.rm = TRUE)
ESS_all$politicisationCGM.x <- ESS_all$politicisation.x - mean(temp_ess$politicisation.x, na.rm = TRUE)

mean(temp_esss$polarization_seats.x, na.rm = TRUE)
ESS_all$polarization_seatCGM.x <- ESS_all$polarization_seats.x - mean(temp_esss$polarization_seats.x, na.rm = TRUE)

mean(temp_ess1$eu_polarisation_ps.x, na.rm = TRUE)
ESS_all$eu_polarisation_psCGM.x <- ESS_all$eu_polarisation_ps.x - mean(temp_ess1$eu_polarisation_ps.x, na.rm = TRUE)

mean(temp_esss2$eu_salience_ps.x, na.rm = TRUE)
ESS_all$eu_salience_psCGM.x <- ESS_all$eu_salience_ps.x - mean(temp_esss2$eu_salience_ps.x, na.rm = TRUE)

mean(temp_essy$politicisation.y, na.rm = TRUE)
ESS_all$politicisationCGM.y <- ESS_all$politicisation.y - mean(temp_essy$politicisation.y, na.rm = TRUE)

mean(temp_ess1y$eu_polarisation_ps.y, na.rm = TRUE)
ESS_all$eu_polarisation_psCGM.y <- ESS_all$eu_polarisation_ps.y - mean(temp_ess1y$eu_polarisation_ps.y, na.rm = TRUE)

mean(temp_esss2y$eu_salience_ps.y, na.rm = TRUE)
ESS_all$eu_salience_psCGM.y <- ESS_all$eu_salience_ps.y - mean(temp_esss2y$eu_salience_ps.y, na.rm = TRUE)

mean(temp_essssss$yearsince, na.rm = TRUE)
ESS_all$yearsinceCGM <- ESS_all$yearsince - mean(temp_essssss$yearsince, na.rm = TRUE)


###################
###MAIN ANALYSIS###
###################

#This code should be adapted according to the different specifications of 
#the model (e.g., remove or add variables, ...)

#LV-2 .X variables refer to CHES Data
#LV-2 .Y variables refer to MARPOR Data

#Remove NA's
ESS_all <- ESS_all%>%drop_na(euftf)
ESS_all$euftf <- as.numeric(ESS_all$euftf)


#MULTILEVEL MODELS WITH CROSS-LEVEL INTERACTIONS
#Economy
model_spl_full_eco <- lmer(euftf ~ # DV
                            stfecoCWC + #LV-1
                            eiscedCWC + #LV-1
                            ageaCWC + #LV-1
                            rural + #LV-1
                            male + #LV-1
                            lrscaleCWC + #LV-1
                            politicisationCGM.x + # LV-2
                            polarization_seatCGM.x + # LV-2
                            IMF + #LV-2
                            Net.payer + #LV-2
                            referendum_EE + #LV-2
                            yearsinceCGM + #LV-2
                            stfecoCWC*politicisationCGM.x + # cross level interaction
                            factor(cntry.x) +
                            factor(year.x) +
                            (1 + stfecoCWC | election_id), # random effect and clustering variables
                          data = ESS_all,
                          na.action = na.omit, 
                          weights = w_pp,
                          control = lmerControl(optimizer = "Nelder_Mead"),  
                          REML = TRUE)  
summary(model_spl_full_eco)
summ(model_spl_full_eco)

#Institutions
model_spl_full_insts <- lmer(euftf ~ # DV
                               trstn_indexCWC + #LV-1
                               eiscedCWC + #LV-1
                               ageaCWC + #LV-1
                               rural + #LV-1
                               male + #LV-1
                               lrscaleCWC + #LV-1
                               politicisationCGM.x + # LV-2
                               polarization_seatCGM.x + # LV-2
                               IMF + #LV-2
                               Net.payer + #LV-2
                               referendum_EE + #LV-2
                               yearsinceCGM + #LV-2
                               trstn_indexCWC*politicisationCGM.x + # cross level interaction
                               factor(cntry.x) +
                               factor(year.x) +
                               (1 + trstn_indexCWC | election_id), # random effect and clustering variables
                             data = ESS_all,
                             na.action = na.omit, 
                             weights = w_pp,
                             control = lmerControl(optimizer = "Nelder_Mead"), 
                             REML = TRUE)  
summary(model_spl_full_insts)
summ(model_spl_full_insts)

#Immigrants
model_spl_full_immig <- lmer(euftf ~ # DV
                                 immig_indexCWC + #LV-1
                                 eiscedCWC +
                                 ageaCWC + #LV-1
                                 rural + #LV-1
                                 male + #LV-1
                                 lrscaleCWC + #LV-1 
                                 politicisationCGM.x + # LV-2
                                 polarization_seatCGM.x + # LV-2
                                 IMF + #LV-2
                                 Net.payer + #LV-2
                                 referendum_EE + #LV-2
                                 yearsinceCGM + #LV-2
                                 immig_indexCWC*politicisationCGM.x + # cross level interaction
                                 factor(cntry.x) +
                                 factor(year.x) +
                                 (1 + immig_indexCWC | election_id), # random effect and clustering variables
                               data = ESS_all,
                               na.action = na.omit, 
                               weights = w_pp,
                               control = lmerControl(optimizer = "Nelder_Mead"), 
                               REML = TRUE)  
summary(model_spl_full_immig)
summ(model_spl_full_immig)


############
###TABLES###
############

setwd()

stargazer(model_spl_full_eco, model_spl_full_insts, model_spl_full_immig, type="text", align= TRUE, 
          covariate.labels = c("Satisfaction National Economy", 
                               "Trust National Institutions", 
                               "Hostility towards Immigrants",
                               "Education", "Age", "Rural Area", "Male",
                               "Left-Right Self-placement",
                               "Party Politicisation EU", 
                               "Party Polarisation (Left-Right)", 
                               "IMF conditionality", 
                               "Net Payer", "Referendum", "Years since membership",
                               "Satisfaction National Economy*Politicisation EU",
                               "Trust National Institutions*Politicisation EU",
                               "Hostility towards Immigrants*Politicisation EU"),
          dep.var.labels = "Attitudes towards European Integration",
          column.sep.width = "1pt",
          model.numbers = FALSE,
          column.labels = c("Model 1", "Model 2", "Model 3", "Model 4", "Model 5"),
          no.space = TRUE, 
          title = "Table 1. Regression Results - MARPOR", 
          single.row = FALSE,
          df = FALSE,
          ci= FALSE,
          out="XXX.tex",
          star.cutoffs = c(0.1, 0.05, 0.001))

############
###GRAPHS###
############

setwd("")

interplot::interplot(model_spl_full_eco, # lmerMod object
                     var1 = "stfecoCWC", # Effect we want to depict
                     var2 = "politicisationCGM.x", # Moderator variable
                     ci = 0.9, # Confidence intervals
                     hist = TRUE) + # Histogram of moderator at bottom
  xlab("Party Politicisation of the EU") +
  ylab("Marginal Effect of Satisfaction with the National Economy on Support towards European Integration") +
  theme_ipsum() +
  geom_hline(yintercept = 0, # Add 0 line of no effect
             linetype = "dashed",
             color = "black",
             size = 1)

ggsave(
  "eco.jpeg",
  plot = last_plot(),
  device = NULL,
  path = NULL,
  scale = 1,
  width = 6,
  height = 7,5,
  units = c("in"),
  dpi = 300,
  limitsize = TRUE
)


interplot::interplot(model_spl_full_insts, # lmerMod object
                     var1 = "trstn_indexCWC", # Effect we want to depict
                     var2 = "politicisationCGM.x", # Moderator variable
                     ci = 0.9, # Confidence intervals
                     hist = TRUE) + # Histogram of moderator at bottom
  xlab("Party Politicisation of the EU") +
  ylab("Marginal Effect of Trust in National Institutions on Support towards European Integration") +
  theme_ipsum() +
  geom_hline(yintercept = 0, # Add 0 line of no effect
             linetype = "dashed",
             color = "black",
             size = 1)

ggsave(
  "insts.jpeg",
  plot = last_plot(),
  device = NULL,
  path = NULL,
  scale = 1,
  width = 6,
  height = 7,5,
  units = c("in", "cm", "mm"),
  dpi = 300,
  limitsize = TRUE
)


interplot::interplot(model_spl_full_immig, # lmerMod object
                     var1 = "immig_indexCWC", # Effect we want to depict
                     var2 = "politicisationCGM.x", # Moderator variable
                     ci = 0.9, # Confidence intervals
                     hist = TRUE) + # Histogram of moderator at bottom
  xlab("Party Politicisation of the EU") +
  ylab("Marginal Effect of Hostility towards Immigrants on Support towards European Integration") +
  theme_ipsum() +
  geom_hline(yintercept = 0, # Add 0 line of no effect
             linetype = "dashed",
             color = "black",
             size = 1)

ggsave(
  "immig.jpeg",
  plot = last_plot(),
  device = NULL,
  path = NULL,
  scale = 1,
  width = 6,
  height = 7,5,
  units = c("in", "cm", "mm"),
  dpi = 300,
  limitsize = TRUE
)

#######################
###Binning Estimator###
#######################


library(base)
ESS_all$election_idd <- gsub('.*^','election',ESS_all$election_id)
ESS_all$w_ppp <- gsub('.*^','a', ESS_all$w_pp)

s1<-cbind.data.frame(stfeco = ESS_all$stfeco, trstn_index = ESS_all$trstn_index, immig_index = ESS_all$immig_index, 
                     euftf = ESS_all$euftf, politicisation.x = ESS_all$politicisation.x,
                     election_id = ESS_all$election_id,
                     eisced = ESS_all$eisced,
                     agea = ESS_all$agea, 
                     rural = ESS_all$rural,
                     male = ESS_all$male,
                     lrscale = ESS_all$lrscale,
                     politicisation.x = ESS_all$politicisation.x,
                     polarization_seats.x = ESS_all$polarization_seats.x,
                     IMF = ESS_all$IMF,
                     Net.payer = ESS_all$Net.payer,
                     referendum_EE = ESS_all$referendum_EE,
                     yearsince = ESS_all$yearsince,
                     cntry.x = ESS_all$cntry.x,
                     year.x = ESS_all$year.x,
                     w_pp = ESS_all$w_pp,
                     w_ppp = ESS_all$w_ppp)

out <- interflex(Y = "euftf", D = "stfeco", X = "politicisation.x", 
                 Z = c("eisced", "agea", "rural", "male",
                       "lrscale", "politicisation.x", "polarization_seats.x",
                       "IMF", "Net.payer", "referendum_EE",
                       "yearsince"), FE = c("cntry.x", "year.x"),
                 data = s1, na.rm = TRUE, estimator = "binning", vcov.type = "cluster", cl = "election_id",
                 weights = c("w_pp"), main = "Marginal Effects", ylim = c(-1, 1))
plot(out,
     theme.bw = TRUE,
     ylab = "Marginal Effect of Satisfaction with the National Economy on Support towards European Integration",
     xlab = "Party Politicisation of the EU",
     cex.lab = 0.7,
     cex.axis = 0.75,
     bin.labs = FALSE,
     CI = TRUE)

out <- interflex(Y = "euftf", D = "trstn_index", X = "politicisation.x", 
                 Z = c("eisced", "agea", "rural", "male",
                       "lrscale", "politicisation.x", "polarization_seats.x",
                       "IMF", "Net.payer", "referendum_EE",
                       "yearsince"), FE = c("cntry.x", "year.x"),
                 data = s1, na.rm = TRUE, estimator = "binning", vcov.type = "cluster", cl = "election_id",
                 weights = c("w_pp"), main = "Marginal Effects", ylim = c(-1, 1))
plot(out,
     theme.bw = TRUE,
     ylab = "Marginal Effect of Trust in National Institutions on Support towards European Integration",
     xlab = "Party Politicisation of the EU",
     cex.lab = 0.7,
     cex.axis = 0.75,
     bin.labs = FALSE,
     CI = TRUE)

out <- interflex(Y = "euftf", D = "immig_index", X = "politicisation.x", 
                 Z = c("eisced", "agea", "rural", "male",
                       "lrscale", "politicisation.x", "polarization_seats.x",
                       "IMF", "Net.payer", "referendum_EE",
                       "yearsince"), FE = c("cntry.x", "year.x"),
                 data = s1, na.rm = TRUE, estimator = "binning", vcov.type = "cluster", cl = "election_id",
                 weights = c("w_pp"), main = "Marginal Effects", ylim = c(-1, 1))
plot(out,
     theme.bw = TRUE,
     ylab = "Marginal Effect of Hostility towars Immigrants on Support towards European Integration",
     xlab = "Party Politicisation of the EU",
     cex.lab = 0.7,
     cex.axis = 0.75,
     bin.labs = FALSE,
     CI = TRUE)

################
###OLS MODELS###
################

model_lm_eco <- lm_robust(euftf ~ stfeco + #LV-1
                            eisced + #LV-1
                            agea + #LV-1
                            rural + #LV-1
                            male + #LV-1
                            lrscale + #LV-1
                            politicisation.x + # LV-2
                            polarization_seats.x + # LV-2
                            IMF + #LV-2
                            Net.payer + #LV-2
                            referendum_EE + #LV-2
                            yearsince + #LV-2
                            stfeco*politicisation.x +
                            factor(cntry.x) + 
                            factor(year.x),
                          data=ESS_all, 
                          weights = ESS_all$w_pp,
                          se = "stata",
                          cluster = election_id)
summary(model_lm_eco)

model_lm_insts <- lm_robust(euftf ~ trstn_index + #LV-1
                              eisced + #LV-1
                              agea + #LV-1
                              rural + #LV-1
                              male + #LV-1
                              lrscale + #LV-1
                              politicisation.x + # LV-2
                              polarization_seats.x + # LV-2
                              IMF + #LV-2
                              Net.payer + #LV-2
                              referendum_EE + #LV-2
                              yearsince + #LV-2
                              trstn_index*politicisation.x +
                              factor(cntry.x) + 
                              factor(year.x),
                            data = ESS_all, 
                            se = "stata",
                            cluster = election_id,
                            weights = ESS_all$w_pp)
summary(model_lm_insts)

model_lm_immig <- lm_robust(euftf ~ immig_index + #LV-1
                              eisced + #LV-1
                              agea + #LV-1
                              rural + #LV-1
                              male + #LV-1
                              lrscale + #LV-1
                              politicisation.x + # LV-2
                              polarization_seats.x + # LV-2
                              IMF + #LV-2
                              Net.payer + #LV-2
                              referendum_EE + #LV-2
                              yearsince + #LV-2
                              immig_index*politicisation.x +
                              factor(cntry.x) + 
                              factor(year.x),
                            data=ESS_all, 
                            weights = ESS_all$w_pp,
                            se_type = "stata",
                            clusters = election_id)
summary(model_lm_immig)

#TABLES
texreg(list(model_lm_eco, model_lm_insts, model_lm_immig), include.ci = FALSE, digits = 3)
screenreg(list(model_lm_eco, model_lm_insts, model_lm_immig), include.ci = FALSE)

#PLOTS
#(It takes some time to run)
margins_model1 <- margins(model_lm_eco, variables = "stfeco", 
                         at = list(politicisation.x = c(7, 10, 20, 30, 40, 50, 80)))
a <- summary(margins_model1, level = 0.90) 

ggplot(data = a, aes(x = politicisation.x, y = AME)) + 
  geom_line() +  
  geom_ribbon(data = a, aes(ymin = lower, ymax = upper), alpha = 0.2) +
  xlab("Party Politicisation of the EU") + 
  ylab("Marginal Effect of Satisfaction with the National Economy on Support towards European Integration") + 
  geom_hline(yintercept = 0, # Add 0 line of no effect
             linetype = "dashed",
             color = "black",
             size = 1) +
  theme_ipsum()

ggsave(
  "ols_eco.jpeg",
  plot = last_plot(),
  device = NULL,
  path = NULL,
  scale = 1,
  width = 6,
  height = 7,5,
  units = c("in", "cm", "mm"),
  dpi = 300,
  limitsize = TRUE
)

margins_model2 <- margins(model_lm_insts, variables = "trstn_index", 
                         at = list(politicisation.x = c(7, 10, 20, 30, 40, 50, 80)))
b <- summary(margins_model2, level = 0.90) 

ggplot(data = b, aes(x = politicisation.x, y = AME)) + 
  geom_line() +  
  geom_ribbon(data = b, aes(ymin = lower, ymax = upper), alpha = 0.2) +
  xlab("Party Politicisation of the EU") + 
  ylab("Marginal Effect of Trust in National Institutions on Support towards European Integration") + 
  geom_hline(yintercept = 0, # Add 0 line of no effect
             linetype = "dashed",
             color = "black",
             size = 1) +
  theme_ipsum()

ggsave(
  "ols_insts.jpeg",
  plot = last_plot(),
  device = NULL,
  path = NULL,
  scale = 1,
  width = 6,
  height = 7,5,
  units = c("in", "cm", "mm"),
  dpi = 300,
  limitsize = TRUE
)




